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PACS 75. 10. Nr - Spin-glass and other random models 

PACS 05 . 70 . Ln - Theories and models of many-electron systems 

PACS 71.55. Jv - Disordered structures; amorphous and glassy solids 

Abstract. - We employ Monte Carlo simulations to investigate the two-time density autocorrela- 
tion function for the two-dimensional Coulomb glass. We find that the nonequilibrium relaxation 
properties of this highly correlated disordered system can be described by a full aging scaling 
ansatz. The scaling exponents are non-universal, and depend on temperature and charge density. 
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Introduction. — Exploring the nature of relaxation 
from out-of-equilibrium initial states towards stationary 
thermal equilibrium has been a major focus of nonequi- 
librium statistical mechanics (for a comprehensive current 
overview, see, e.g., Ref. pQ). Fundamental questions con- 
cern the scaling properties of two-time correlation func- 
tions in the so-called aging regime (to be discussed more 
extensively below), namely: (1) What are the basic scal- 
ing forms to describe numerical and experimental data? 
(2) Are the ensuing scaling exponents and scaling func- 
tions universal? (3) Which physical properties determine 
the scaling features in the aging regime? Aside from con- 
tributing to our still rather incomplete understanding of 
processes far from equilibrium, providing comprehensive 
and reliable answers to questions (2) and (3) could have 
important ramifications to materials science as well. For 
example, a clear picture of which features of nonequilib- 
rium relaxation processes contain specific information on 
the system under investigation in contrast to which prop- 
erties are generic may provide a novel method of sample 
characterization. To this end, various model systems have 
been carefully investigated through numerical simulations; 
a few simple models even allow for analytical treatments 



Disordered semiconductors in the vicinity of a metal- 
insulator transition constitute an intriguing experimental 
system displaying complex relaxation phenomena. Re- 
laxation measurements for the conductivity of a two- 
dimensional silicon sample have yielded unambiguous ev- 
idence for aging effects [2], see also Ref. [5]. The interpre- 
tation of these experiments continues to build upon the 



early theoretical predictions of Coulomb glass models [3HS] 
that incorporate the essential physical aspects of charge 
carriers localized at random positions in space and inter- 
acting through long-range forces. Research into the non- 
trivial equilibrium properties [71415] and non-equilibrium 
relaxation phenomena [16h24] of the Coulomb glass over 
the past two decades have considerably advanced our un- 
derstanding of this paradigmatic model system for highly 
correlated disordered materials. Our goal in this paper 
is to examine the fundamental questions (l)-(3) for the 
Coulomb glass in two dimensions, with the aim to better 
understand universality aspects and scaling properties of 
disordered semiconductors. (We note that a very similar 
model, with the Coulomb potential essentially replaced 
by a logarithmic repulsion, describes the Bose glass phase 
of magnetic flux lines in type-II superconductors that are 
pinned to columnar defects, see Refs. [2"5h27j .) Specifi- 
cally, wc address the detailed scaling form for the two-time 
density autocorrelation function, and the dependence on 
temperature and carrier density of the non-equilibrium re- 
laxation kinetics. Our work builds upon and extends the 
investigations of Grcmpel et al. who employed a very sim- 
ilar Monte Carlo simulation method [2UJ|2T]; however, we 
decided to extend our studies to multiple carrier densi- 
ties in addition to varying temperatures. This allows us 
to study in a very systematic way the aging properties 
of the two-dimensional Coulomb glass, utilizing different 
scaling forms for the density autocorrelation function (see 
also Ref. [35J in the context of a spin glass model). 

The Coulomb glass model and Monte Carlo sim- 
ulation procedure. — Our basic model system is a 
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Coulomb glass in two dimensions. The Coulomb glass 
model, introduced by Efros and Shklovskii pjj, consists 
of multiple localized pinning sites available to the charge 
carriers. Because of the strong intra-site correlations these 
sites can only contain a single charge carrier at most. The 
system is dominated by long-range repulsive Coulomb in- 
teractions V(r) ~ 1/r and the spatial disorder is induced 
by the randomly located (on a continuum) available sites. 
The Hamiltonian of the Coulomb glass model reads [4]-[6] 



H 



E 



Tliipi 
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R.i — Ri 



(1) 



where R^ , ifi, and n, respectively denote the position vec- 
tor (here in two dimensions), (bare) site energy, and occu- 
pancy of the ith site, i = 1, . . . , N. The occupancy m can 
only take on the values or 1. The first term corresponds 
to (random) site energies assigned to each accessible lo- 
cation; since the system is dominated by the long-range 
forces, we chose ipi = to further simplify the model, while 
still allowing the positions Ri to be continuous and ran- 
dom [9][T7j[20]. (Alternatively, one could allow ifi to take 
on random values while the positions Ri would be discrete 
and set on a lattice.) The second contribution encapsu- 
lates the repulsive Coulomb interactions (with dielectric 
constant k). In order to maintain global charge neutral- 
ity, a uniform relative charge density K = J2i n i/N is 
inserted; K can also be described as the total carrier den- 
sity per site, or filling fraction. We note that with (fi = 
the Hamiltonian ([1} displays particle-hole symmetry, i.e., 
systems with K = 0.5 + k and K = 0.5 — k are equivalent. 
Upon replacing the occupation numbers with spin vari- 
ables <Tj = 2rii — 1 = ±1, one sees that the Coulomb glass 
maps onto a random-site random-field antifcrromagnctic 
Ising model with long-range exchange interactions [TJ. 

We employ a Monte Carlo simulation algorithm to at 
least heuristically model the kinetics in the Coulomb glass 
20.21]. At each time step, one randomly selected charge 
carrier attempts to hop from an occupied site a to an 
empty site b. Two factors determine the hopping success 
rate, namely a strongly distance-dependent tunneling pro- 
cess (in semiconductors mediated through phonons) and 
thermally activated jumps over energy barriers |20j . i.e., 



r a ^ & = r - 1 e - 2 ^^min[l, 



,- AE ab / T] ; (2) 



where tq represents a microscopic time scale, R a b = 
R a — R(,| is the distance between sites a and o, £ charac- 
terizes the spatial extension of the localized carrier wave 
functions, and we have set ks = 1. The first, spatially ex- 
ponential term in ([2]) is derived from quantum tunneling, 
while the second exponential term is due to thermodynam- 
ics: A thermally activated hop from sites a to b entails the 
energy difference 



with the (interacting) site energies 

n b - K 
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In the following, distances are measured relative to the 
average separation between sites ao, and energies as well 
as temperature scales are measured relative to the typical 
Coulomb energy Ec — e 2 / kclo [20]. We set the spatial 
extension of the localized wave functions, £, to ao, as in 
Refs. [2U1I2TJ . We have explored other values for £ as well, 
and (within the applicability range of the model) found 
the ensuing results to simply scale with To; hence £ = ao 
was chosen for simplification. 

The Monte Carlo simulations were initiated by ran- 
domly placing N sites within a square simulation cell of 
length L containing N e — KN charge carriers, where 
N = L 2 . We performed simulations for systems with 
L = 8, 10, 16, 32; with temperatures in the range of 
0.01 < T < 0.05; and carrier densities in the range of 
0.375 < K < 0.5 (also equivalent to 0.5 < K < 0.625 due 
to particle-hole symmetry). Running the simulations with 
various system sizes L, we noticed no measurable finite- 
size effects, as demonstrated in Fig. 2(a) Temperatures 
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larger than 0.03 turned out not to be useful for our study 
of aging processes since equilibrium was then reached far 
too quickly. In contrast, for T < 0.01, the kinetics slowed 
down too much for gathering statistically significant data 
within computationally reasonable time frames. As will be 
discussed in more detail below, the dynamics also freezes 
out within the numerically accessible simulation times for 
K < 0.4 (or K > 0.6). 

Periodic boundary conditions were used and the poten- 
tial due to charges outside the cell was calculated by mir- 
roring the simulation cell on each side. Initially, the charge 
carriers were placed randomly at available sites to mimic 
a quench from very high temperatures. Then the system 
was evolved at temperature T with the dynamics described 
by the generalized Metropolis rate ^. The thermally ac- 
tivated, tunneling-controlled (variable-range) hopping al- 
gorithm begins with randomly choosing an occupied site 
a and assigning a normalized probability proportional to 
exp(— 2R a b/t;) to each empty site. An empty site b is then 
chosen from this probability distribution and a hop from 
a to b is attempted. The success probability for this move 
is determined by min[l, e~ AEab / T }. There are iV hop at- 
tempts for each Monte Carlo step (MCS). Our simulation 
runs consisted of 1 x 10 6 MCS and results were averaged 
over 3000 different realizations of the initial conditions and 
of the random number sequences. 

In thermal equilibrium, the long-range correlations and 
spatial disorder conspire to produce a soft gap in the (in- 
teracting) density of states g(e), as famously first noted by 
Efros and Shklovskii [J] . The repulsive interactions induce 
strong spatial anticorrelations, which in turn suppress the 
availability of states with energies close to the chemical 
potential /j c that at T = separates the occupied and 
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Fig. 1: Coulomb gap formation in a Coulomb glass (L = 16, 
iV = 256 sites) at half filling K = 0.5 and at temperature 
T = 0.02. The plot shows how g(fJ, c ) quickly approaches zero. 
The density of states is displayed in the inset at t — 1 MCS 
(dashed) and t = 100 MCS (solid), when the Coulomb gap 
correlations are clearly established. 



empty energy levels. The ensuing Coulomb gap appears to 
asymptotically be described by a power law g(e) ~ e— fi c \ s 
[5UT4]: within a simple mean- field approach, s = (d/a) — 1 
in d dimensions for a system with long-range repulsive po- 
tential V(r) ~ l/r°\ We have monitored the emergence of 
the Coulomb gap in our Monte Carlo simulations. Quite 
remarkably, the soft gap forms very fast, as illustrated 
in Fig. [TJ which shows how the density of states at the 
chemical potential g{fi c ) quickly approaches zero. Indeed, 
the Coulomb gap is already pronounced after as few as 
10 MCS, and clearly established in less than 100 MCS. 
In the following, we proceed to carefully analyze the slow 
decay of local density correlations; we remark that the in- 
teresting aging kinetics happens in a time period when the 
Coulomb gap itself is already well-established and there- 
fore addresses slow dynamics in a highly correlated disor- 
dered system. 

Scaling properties or the two-time density auto- 
correlation function in the aging regime. In the 

following, we focus on the (normalized) two-time carrier 
density autocorrelation function |20j , 



C(t,s) 



(ni(t)ni(s)) - 
K{l-K) 



J2ini(t)ni(s)-NK 
NK(1 - K) 



-2 



(5) 

where t > s, and the waiting time s refers to the time 
elapsed since the high-temperature quench (initiation of 
the Monte Carlo simulations). If both s and t are large 
compared to microscopic time scales tq , and also well sep- 
arated, i.e., if 



t > T , s 3> r , t - s > t , 



(6) 



timc-translational invariance is broken and C depends on 
both t and s separately. In the aging scaling regime, one 
then expects the autocorrelation function to obey the gen- 



eral scaling form (we follow the notation in Ref. Q~]), 

C(t,s) = s- b f c (t/ S n- (7) 

For certain simple systems obeying purely relaxational dy- 
namics, namely the one-dimensional kinetic Ising model 
[5^1130], the time-dependent Ginzburg-Landau models 
quenched to a critical point, [22 [32 and the spherical 
model A coarsening dynamics in the low-temperature 
phase, with short-range [331134] or long-range [35H37J inter- 
actions, the aging scaling regime is amenable to analytic 
treatment. In these cases one finds that the ansatz ([7]) is 
satisfied with exponent fj, = 1, a situation commonly re- 
ferred to as full aging scaling. For more complex systems 
such as spin glasses the possibility of a subaging scaling 
with n < 1 [1] has been discussed in the literature. For 
fi = 1 and t/s — > oo the full aging scaling function follows 
a power law [T], 

i-Ac/2 



fc(t/s) ~ (t/s)- 



(8) 



with the autocorrelation exponent Ac* and the dynamic 
exponent z. Overall, this presents us with three scaling 
exponents: b, fi, and Xc/z. 

From the averaged charge density autocorrelation func- 
tion (0 we extracted the scaling exponents using an inter- 
polation method motivated by Ref. [38]. The total vari- 
ance, or spread, of the data was calculated by comparing 
the data points to polynomial interpolations of all other 
curves. This variance was only measured within the scal- 
ing regime. We have chosen to either assume full aging 
scaling, fi = 1 in Eq. ([7]), and applied the polynomial in- 
terpolation method to obtain the scaling exponent b > 0; 
or to set 6 = (simple scaling ansatz) and alternatively 
determine the subaging exponent fi < 1. 

Sample sets of results (with L = 16, N = 256 sites) 
are shown in Figs. [5] and [31 The density autocorrelation 
data at half filling K = 0.5 displayed in Fig. 2(a) match 
nicely with the results of Ref. [20]; however, the scaling 



forms used in Figs. 2(b) and 2(c) are different from those 

~ the 



2(a) 



employed by Grcmpel. As is apparent from Fig. 
density autocorrelation function is characterized by two 
temporal regimes. In the first time regime, the density au- 
tocorrelation relaxes towards a plateau (often referred to 
as j3 relaxation in the glass literature) , whereas in the sec- 
ond time regime the system very slowly approaches equi- 
librium (a relaxation) [SHIED]- This later regime displays 
breaking of time translation invariance and aging scaling, 
provided the inequalities ([6]) hold. With our current algo- 
rithm and available hardware, we cannot effectively inves- 
tigate systems with densities lower than K = 0.40625 and 
temperatures lower than T = 0.02, since in these configu- 
rations the plateau regions persist for much longer times 
than arc computationally accessible for us, i.e., the system 
essentially freezes into a structure that corresponds to the 
/3 relaxation plateau. Notice that the autocorrelation data 
obtained for a smaller system (L = 10, N = 100 sites) co- 
incide with those for the larger system, demonstrating the 
absence of finite-size effects in the observed time window. 
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Fig. 2: Carrier density autocorrelation function C(t, s) at half 
filling K = 0.5, in a system of linear extension L = f 6 (N = 256 
sites) at temperature T = 0.02, measured for various waiting 
times s = 10 2 ,2 x 10 2 ,5 x 10 2 ,10 3 ,2 x 10 3 ,5 x 10 3 ,10 4 (from 
bottom to top), (a) The plot vs. t — s demonstrates the break- 
ing of time-translation invariance. The • symbols represent 
data obtained from a system of length L = 10 (N — 100 sites), 
(b) A full aging scaling plot, which assumes fj, — 1, yields the 
scaling exponents b — 0.038 and Xc/z — 0.103. (c) A subaging 
scaling analysis, fixing 6 = 0, gives fi = 0.67. 



(c) 



Fig. 3: (a) Carrier density autocorrelation function C(t, s) at 
carrier density K = 0.4375 (again with L — 16, N — 256, and 
T = 0.02), measured for various waiting times s = 10 2 ,2 x 
10 2 ,5 x 10 2 ,10 3 ,2 x 10 3 ,5 x 10 3 ,10 4 (from bottom to top). 
Note the markedly slower decay as compared to the date at half 
filling shown in Fig. [5] (b) A full aging scaling plot (^ = 1) 
gives b = 0.006 and \c/z — 0.075. (c) A subaging scaling 
analysis (b = 0) results in /i = 0.90. 
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Table 1: Charge carrier density dependence of the full aging 
scaling exponents b and Xc/z (with /x = 1 and T — 0.02). 



K 



Xc/z 



0.40625 0.001 ± 0.002 0.046 ± 0.004 

0.4375 0.006 ±0.002 0.075 ± 0.004 

0.46875 0.021 ±0.004 0.089 ±0.003 

0.5 0.038 ±0.006 0.103 ±0.004 



ponents b and Xc/z (with /j, — 1 and K — 0.5). 

T b Xc/z 



0.01 
0.02 
0.03 



-0.001 ±0.001 0.036 ±0.004 
0.038 ±0.006 0.103 ±0.004 
0.080 ±0.008 0.173 ±0.005 



Table 3: Charge carrier density dependence of the subaging 
exponent fj, (with b = and T — 0.02). 



K 



/' 



0.40625 




1.00 ±0.01 


0.4375 




0.90 ±0.01 


0.46875 




0.78 ±0.01 


0.5 




0.67 ±0.02 


Table 4: Temperature 


dependence 


of the subaging 


(with b = and K = 


.5). 




T 




A* 


0.01 




0.98 ±0.06 


0.02 




0.67 ±0.02 


0.03 




0.54 ±0.02 



In Fig. |2(b)| we plot the same data as in Fig. 2(a) af- 
ter applying a rescaling procedure that assumes full aging 
(where /i = 1). One notes that the systematic deviations 
for small t become less and less important for increasing 
waiting times s. From the scaling ansatz ([7]), which is 
well founded theoretically [I], the scaling exponent b can 
be reliably extracted from the data sets corresponding to 
larger waiting times. As expected from Eq. ((5|), a power 
law is observed at long times, which in turn yields the 
scaling exponent Xc/z. At half filling K = 0.5, we find 
b = 0.038 ± 0.006 and X c /z = 0.103 ± 0.004. Alterna- 
tively, we may employ a subaging scaling ansatz, setting 
6 = 0. Figure 2(c) shows the results from scaling using fi 



only. Also in this case deviations are observed for t small, 
but an improvement is less apparent when increasing the 
waiting time. Collapsing the data on a single master curve 
at large scaling arguments works apparently about equally 
well with either method employed in (b) or (c). In fact, al- 
lowing nontrivial values for both scaling exponents, 6 > 
and /i < 1, yields comparable scaling collapse quality. 
However, given that theoretical analysis (admittedly for 
comparatively simple systems) yields the full aging sce- 
nario with fi — 1 and generally b > 0, we would argue 
that subaging scaling, as used in previous studies [2TJ1I21] . 
need not be invoked to describe the nonequilibrium relax- 
ation phenomena in the Coulomb glass. 

Figure [3] shows the results when the same scaling meth- 
ods are applied to our data for relative charge carrier 
density K = 0.4375 (or K = 0.5625). As the carrier 
density deviates from half filling, the relaxation becomes 
markedly slower. The full aging scaling ansatz (/x = 1) 
yields an excellent data collapse for all s and t values. 
The same conclusion holds for other charge carrier den- 
sities K 7^ 0.5. For K = 0.4375, from the scaling anal- 
ysis based on Fig. |3(b)| we infer the rather low value 
b = 0.006 ± 0.002 along with X c / z = 0.075 ± 0.004. If we 
instead assume b = 0, the subaging scaling from Fig. 3(c) 
gives /i = 0.90 ± 0.01. Clearly, the aging scaling expo- 



nents arc not universal in the Coulomb glass, but strongly 
depend on the charge carrier density. At least near half- 
filling, in the time window that is accessible to us we ob- 
serve power-law scaling, albeit with numerically small ex- 
ponent values b and Xc/z, different from, but close to the 
logarithmic dependence obtained by the mean-field type 
analysis in Refs. [2"2"ll2"5] . 

We ran sets of 3000 simulation runs each at various val- 
ues of T and K to systematically explore the dependence 
of the aging scaling exponents on the temperature (already 
noted in Refs. (2HH2]) and carrier density. The resulting 
exponent values for b and Xc/z as obtained within the 
full aging framework (setting /j, = 1) are listed in Tables Q] 
and [51 Physically, smaller values of the exponents imply 
slower relaxation kinetics. The trend that is observed is 
that as density drops away from half filling, the relaxation 
slows. A similar trend is seen in the temperature depen- 
dence; namely, the relaxation processes become slower at 
lower temperatures, as one would naturally expect. If, on 
the other hand, we enforce 6 = and instead allow for 
the subaging scaling scenario, we find the exponent values 
fj, listed in Tables [3] and |H Note that a larger value of 
/z leads to slower relaxation. The dependence of the sub- 
aging exponent thus follows the same qualitative trends 
as observed in full aging scaling. Naturally, as in the full 
aging scaling analysis the exponent b approaches zero at 
low temperatures or away from half-filling, in the subag- 
ing scaling analysis /j, — > 1. (In fact, for cither T = 0.01 or 
K = 0.40625 our simulations approach the aforementioned 
limitations in run times; we list the corresponding expo- 
nent values with large relative errors nevertheless, since 
they still confirm the general trends.) We remark that our 
findings for the scaling of nonequilibrium relaxation pro- 
cesses in the Coulomb glass align well with a recent study 
of the aging kinetics in the two-dimensional ferromagnetic 
Ising model with uniform bond disorder: Ref. |41j also 
finds non-universal scaling exponents, depending on the 
ratio of disorder distribution width and temperature. 
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Conclusion. — We have investigated nonequilibrium 
relaxation processes of the two-dimensional Coulomb glass 
model at low temperatures via Monte Carlo simulations, 
and confirmed that its scaling features are governed by 
the simple aging scaling form. Through the use of polyno- 
mial interpolation methods, three scaling exponents were 
extracted and found to follow a common trend: as either 
the temperature decreases or the charge carrier density de- 
viates more from half-filling, the exponents reflect slower 
relaxation kinetics. The aging scaling exponents arc thus 
found to be non-universal, as is the case for disordered 
magnetic systems [4"I] . Our results indicate that the in- 
clusion of a subaging exponent /i ^ 1 is unnecessary since a 
nonzero scaling exponent b encapsulates the same charac- 
teristics, and the corresponding simple aging scaling form 
is more firmly grounded on theoretical analysis. 

We are presently working on extending the scope of our 
study into three dimensions. Preliminary data for the den- 
sity autocorrelation function scaling in the aging regime 
suggest the same basic features and trends that we observe 
in two dimensions [55] • We also plan to implement loga- 
rithmic repulsion instead of the 1/r potential used in the 
Coulomb glass, with the aim to study aging phenomena 
in the Bose glass phase of type-II superconductors with 
columnar defects [25H27] , Our ultimate goal will be to 
more closely relate our observables and findings to exper- 
imental setups and measurements. 
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